function [newSeq newSeqString] = addZerosAndTargets(baseSeq, SOA, numTargets, blankBlocks)

newSeqString = blanks(length(baseSeq));
newSeqString(newSeqString==' ') = '-';

numLabels = 2;
seqLength = length(baseSeq);
wordLength = log2(seqLength)/log2(numLabels); % logB(x) = log2(x)/log2(B), in case numLabels ~= 2

TRDuration = 3000;

numSyncTRs = lcm(SOA, TRDuration) / TRDuration; % Number of TRs required to sync scanner-experiment
numSyncTrials = lcm(SOA, TRDuration) / SOA % Number of trials required to sync scanner-experiment
numTrials = ceil(seqLength/numSyncTrials)*numSyncTrials;
numTRs = ceil(seqLength/numSyncTrials)*numSyncTrials*SOA/TRDuration;


%% Wrap sequence around so that the number of trials is a multiple of numSyncTrials
trialsToComplete = numTrials-seqLength; % Number of trials to be taken from the end of the sequence and include in the beginning

for i=1:trialsToComplete
    newSeqString = ['B' newSeqString];
end


%% Re-organize sequence in groups of numSyncTrials
newSeqString = reshape(newSeqString,numSyncTrials,numTrials/numSyncTrials);


%% Add target trials

targetPositions = sort(ceil(rand(1,numTargets)*(numTrials/numSyncTrials - 1))); % Target blocks can be added between any two blocks, but not in the beginning/end of the sequence

targetBlock = [];
for i=1:numSyncTrials
    if i==1
        targetBlock = [targetBlock;'T'];
    else
        targetBlock = [targetBlock;'R'];
    end
end

for i=1:length(targetPositions)
    lineIndex = targetPositions(i) + 1 + (size(newSeqString,2)-numTrials/numSyncTrials);
    newSeqString = [newSeqString(:,1:(lineIndex-1)) targetBlock newSeqString(:,lineIndex:end)];
end


%% Add blank blocks

if blankBlocks > 1
    blankPositions = floor(0:(size(newSeqString,2)/(blankBlocks-1)):size(newSeqString,2));
else
    blankPositions = 0;
end

blankBlock = [];
for i=1:numSyncTrials
    blankBlock = [blankBlock;'B'];
end

for i=1:length(blankPositions)
    lineIndex = blankPositions(i) + 1 + (size(newSeqString,2)-(numTrials/numSyncTrials + numTargets));
    if lineIndex <= size(newSeqString,2)
        newSeqString = [newSeqString(:,1:(lineIndex-1)) blankBlock newSeqString(:,lineIndex:end)];
    else
        newSeqString = [newSeqString(:,1:(lineIndex-1)) blankBlock];
    end
end


%% Set regions of overlap

overlapBlock = blanks(wordLength-1);
overlapBlock(overlapBlock==' ') = 'O';

for i = 1:length(newSeqString(:))
    if newSeqString(i)=='-' && newSeqString(i-1)~='-'
        newSeqString(i-(wordLength-1):i-1) = overlapBlock;
    end
end


%% Now, create the real sequence
newSeq = nan(size(newSeqString));

trialIndex = 1;

% Add valid trials
for i=1:length(newSeq(:))
    if newSeqString(i) == '-'
        newSeq(i) = baseSeq(trialIndex);
        trialIndex = trialIndex + 1;
    end
end

% Add overlap trials
copyingOverlapStarted = 0;
j = nan;
for i=1:length(newSeq(:))
    if ((newSeqString(i) == 'O') && (copyingOverlapStarted == 0))
        copyingOverlapStarted = 1;
        j=i;
        % Find the last valid trial
        while newSeqString(j) ~= '-'
            if j > 1
                j = j - 1;
            else
                j = length(newSeq(:));
            end 
        end
        j = j - (wordLength-2);
        newSeq(i) = newSeq(j);
        j = j + 1;
    elseif ((newSeqString(i) == 'O') && (copyingOverlapStarted == 1))
        newSeq(i) = newSeq(j);
        j = j + 1;
    else
        copyingOverlapStarted = 0;
        j = nan;
    end
end

for i=1:length(newSeq(:))
    switch newSeqString(i)
        case 'B'
            newSeq(i) = 0;
        case 'T'
            newSeq(i) = -1;
        case 'R'
            newSeq(i) = ceil(rand*numLabels); % Random number between 1 and 2
    end
end
    
    































